
addpath ../Generator
addpath ../Generator/GrTheory
addpath ../PDCO
addpath ../CRNTSolvers

%Generate a random single linkage class and try 
%to converge to the fixed point of the mapping
if(~isvarname('Y'))
		m = 9;
		n = 15;
		r = 3;
		l = 2;
		%Set the seed to be able to reproduce the experiment
		RandStream.setDefaultStream(RandStream('mt19937ar','seed',0));
		%Generate a random set of n complexes on n species, each with at most 3 species
		Y = YGenerator(m,n,r);
		%Generate the strongly connected graph with a single linkage class`
		Ak = AkGenerator(n,0.2,l);
else
		m = size(Y,1);
		n = size(Y,2);
end
A = Ak' - diag(diag(Ak));
D = -diag(diag(Ak));

YAt = Y*A';
YD = Y*D;

gamma0 = sum(Y*diag(D));
gamma_vec = gamma0*(0.2:0.2:5);
C_vecs = [];
for gamma = 0.2:0.2:5
    x0 = gamma*ones(n,1);
    [iter,v_vecs,lmda_vecs,mass_infeas,mass_action_infeas]=...
        SolverFpIterationHomPDCO(Y,Ak,1e-6,500,0.5,x0);
    C_vecs  = [C_vecs exp(lmda_vecs(1:m,end))];
end

figure(1)
imagesc(sign(Ak))

figure(2)
title('Equilibrium Concentrations vs \gamma_{0}');
ylabel('Equilibrium Concentrations of species');
xlabel('\gamma_{0}');
for i = 1:m
    subplot(3,2,i), plot(gamma_vec, C_vecs(i,:));
end

C = exp(lmda_vecs);
figure(3)
title('Concentration paths of species vs iteration of fixed point algorithm');
ylabel('Concentration of species');
xlabel('Iteration');
for i = 1:m
    subplot(3,2,i), plot(C(i,:));
end
